#!/bin/csh
# 2d_plots.csh
#############################################################################
# 2D contour plots of YACC output fields
#
# - Set the simulation directory containing YACC's output 
#   in the file "simulation.dir"
# - To make a dimensional plot, prepare the file scaling.dat in the simulation directory
# - Set the options in this file (see below)
# - Run the script using:
#    ./2d_plots.csh n field xmin xmax zmin zmax    
#    - 'n' is the time step to be plotted 
#      Note that n must be given in the following form:
#      000001 or 000010 or 000100 or 001000 etc.
#
#    - 'field' is the field to plot. Choose among:
#      temperature
#      tempanomaly: temperature anomaly with respect to the laterally averaged profile
#      tbuoyancy: thermal buoyancy (as tempanomaly but including expansivity
#      velocity: RMS velocity
#      viscosity
#      viscdiss: viscous dissipation 
#      composition
#      companomaly: composition anomaly with respect to the laterally averaged profile 
#      conductivity: thermal conductivity
#      expansivity: thermal expansicity
#      heats: Rayleigh number for internal heating
#      depletion
#      solidus
#      strain: second invariant of the strain rate tensor
#      defmech: deformation mechanism
#     
#     - 'xmin', 'xmax', 'zmin', 'zmax' (optional) are the coordinates of the 
#       domain limits within which the contour plot is needed. If dimensional_plot = 1, the coordinates 
#       must be expressed in km. If nothing is provided the whole domain is plotted.
#########################################################################################

#--------------------------------------------------------------------------------------
# Options
#--------------------------------------------------------------------------------------
# dimensional_plot: if 1 the file scaling.dat must be in the simulation directory, otherwise 
#                   a non-dimensional plot is made (default: 0)
# makejpg: make also jpg files (default: 1)
# jpgdensity: resolution of the converted jpg files (default: 200)
# deletepdf: delete pdf files (done only if makejpg=1, default: 1)
# dimz: height of the plot in cm (default: 5 cm)
# velocity_vectors: plot velocity vectors (default: 1)
# const_vec_leng: 1 length of the vectors is constant, 
#                 0 length of the vectors proprotional to the velocity field
#                 (default: 1)
# vec_col: RGB color of the velocity vectors (default: 50/50/50, i.e. gray)
# deltaux, deltauz: non-dimensional intervals along x an z for plotting velocity vectors (default 0.05)
# streamlines: plot streamlines (default: 0)
# stream_col: RGB color of the streamlines (default: 100/100/100, i.e. light gray)
# deltat, nsource, nsteps: streamlines parameters (time-interval for particle advection,
#   number of source points, number of time steps)
# plot_mesh: overplots the underlying mesh (default: 0)
# write_time: writes the current simulation time as title (default: 1)
# phase1, phase2, ...: plot the corresponding phase boundary (default: 0)
# phase_col: RGB color of the phase boundaries (default: 50/50/50, i.e. gray)
# X0: x-coordinate in cm of the bottom-left corner of the plot (default: 3 cm)
# Y0: Y-coordinate in cm of the bottom-left corner of the plot (default: 5 cm)
# xi, zi: non-dimensional intervals along x and z for the contour plot (default: 0.002)
# dirplots: directory where to store the plots (default: $dirout/Plots/$field)
#--------------------------------------------------------------------------------------
set dimensional_plot = 0
set makejpg = 1
set jpgdensity = 200
set deletepdf = 1
set dimz = 5
set velocity_vectors = 1
set const_vec_leng = 0
set vec_col = `echo "50/50/50"`
set deltaux = 0.05
set deltauz = 0.05 
set streamlines = 0
set stream_col = `echo "100/100/100"`
set deltat = 1d-5
set nsource = 100
set nsteps = 100
set plot_mesh = 0
set write_time = 1
set phase1 = 0
set phase2 = 0
set phase3 = 0
set phase4 = 0
set phase5 = 0
set phase_col = `echo "250/250/250"`
set X0 = 3
set Y0 = 5
set xi = 0.002
set zi = 0.002
#--------------------------------------------------------------------------------------

#-----------------------------------------------
#  Read simulation directory and scaling file
#-----------------------------------------------
source simulation.dir

if (! -d $dirout) then
  echo $dirout 'does not exist. Check the path stored in the variable "dirout" '
  exit
endif

#-------------------------------------------------------------------
# If the scaling file is not found, non-dimensional plots are made
#-------------------------------------------------------------------
if ($dimensional_plot == 1) then
  if (! -f {$dirout}scaling.dat) then
    echo {$dirout}scaling.dat 'for dimensional plot not found.'
    set dimensional_plot = 0
  endif
endif
  
if ($dimensional_plot == 1) then
  echo 'Dimensional plot'  
  source {$dirout}scaling.dat
else
  echo 'Non-dimensional plot'  
  set D_scale = 1.0
  set T_scale = 1.0
  set T_surf = 0.0
  set eta_ref = 1.0
  set k_ref = 1.0
  set alpha_ref = 1.0
  set time_scale = 1.0
  set V_scale = 1.0
  set eps_scale = 1.0
  set heats_scale = 1.0
endif

set deltaux = `echo "$deltaux * $D_scale" | bc -l`
set deltauz = `echo "$deltauz * $D_scale" | bc -l`
set xi = `echo "$xi * $D_scale" | bc -l`
set zi = `echo "$zi * $D_scale" | bc -l`

#--------------------------------------------------------------------------------------
# Specific options are set separately for each field (see below):
#--------------------------------------------------------------------------------------
# xlabel: label for x-axis (default: none) 
# zlabel: label for z-axis (default: none)
# clabel: label for color scale (default: field)
# cscale: gmt color scale (polar, seis, wysiwyg, rainbow, jet, etc.)
# ctick: tick intervals for color scale
# contour_lines: contour lines over color plot
# stepcontour: interval for contour lines
# opt: additional options for gmt color scale (
# adaptive_cscale: 1 set color scale range automatically
#                  0 set color scale range manually according to min max 
# min: set manually the minimum of the field 
# max: set manually the maximum of the field
# pdfoutput: name of pdf output file (default field_cplot_{$tstep}.pdf)
# jpgoutput: name of jpg output file (default field_cplot_{$tstep}.jpg)
#--------------------------------------------------------------------------------------

set dirplots = {$dirout}Plots/ 

#-------------------------------
# Select YACC's output files
#-------------------------------
if (! -d $dirplots) then
  mkdir $dirplots
endif

set tstep = ${1} 
set nchar = `echo $tstep | awk '{print length($0)}'`
if ($nchar != 6) then
  echo "Wrong format of the input time step."
  exit
endif

set gridfile = ${dirout}grd_{$tstep}_.dat
set profile = ${dirout}profiles_{$tstep}_.dat
if (! -f $gridfile || ! -f $profile ) then
  echo 'Either ' $gridfile ' or ' $profile ' or both do not exist.'
  exit
endif
echo "YACC grid output file:" $gridfile
echo "YACC profile output file:" $profile

if($plot_mesh == 1) then
  set meshfile = ${dirout}mesh.dat
  if (! -f $meshfile) then
    echo $meshfile ' does not exist. The mesh will not be plotted.'
    set plot_mesh = 0
  endif
endif

if($streamlines == 1) then
  set streamfile = streamlines_{$tstep}_.dat
endif

#------------------------------------------------------------
# Read and set current simulation time in the variable ctime
# (it can be used e.g. as plot title)
#------------------------------------------------------------
set ctime = `awk 'NR == 1 {print $5}' $gridfile`
if($dimensional_plot == 1) then
  set ctime = `awk 'BEGIN {print "'"$ctime"'"*"'"$time_scale"'"}'`
  set ctime = `echo $ctime Ma`
endif

#--------------------------------------------------------------------------------------------
# Read and set the aspect ratio 
# NB: it will NOT work with old input files which did not contain aspect ratio information. 
#     in that case it must be set manually 'set ar = ...'
#--------------------------------------------------------------------------------------------
set ar = `awk 'NR == 1 {print $11}' $gridfile`
if ($ar == "") then
 echo 'You are using an old output file with no information on the aspect ratio.'
 echo 'Set the aspect ratio explicitly in 2d_plots.csh as follows: "set ar = ... "'
 exit
endif

#------------------------------------------------------
# Read the field to plot from standard input
#------------------------------------------------------
set field = ${2}
if (${2} !~ temperature && ${2} !~ tempanomaly && ${2} !~ tbuoyancy && ${2} !~ velocity && ${2} !~ viscosity && \
    ${2} !~ viscdiss && ${2} !~ composition && ${2} !~ companomaly && ${2} !~ conductivity && ${2} !~ expansivity \
    && ${2} !~ heats && ${2} !~ depletion && ${2} !~ solidus && ${2} !~ strain && ${2} !~ defmech) then
  echo 'The field you want to plot is either not (yet) supported in this routine or was misspelled'
  exit
endif

if (! -d $dirplots$field) then
  mkdir $dirplots$field
endif

#------------------------------------------------
# Set horizontal and vertical range to plot
# default is whole domain, i.e. x=[0,ar] z=[0,1]
#------------------------------------------------
if (${3} != "" && ${4} == "" || ${4} != "" && ${5} == "" || ${5} != "" && ${6} == "") then
  echo 'To plot only part of the domain, xmin, xmax, zmin and zmax must all be set.'
  exit
endif
if (${3} == "") then
  set xmin = 0
  set xmax = `echo "$ar * $D_scale" | bc -l`
  set zmin = 0
  set zmax = `echo "$D_scale" | bc -l`
else
  set xmin = ${3} 
  set xmax = ${4}
  set zmin = ${5}
  set zmax = ${6}
endif  

#----------------------------------------------------
# Set size of the plot and of the color bars (in cm)
#----------------------------------------------------
set dimx = `echo "($xmax-$xmin)/($zmax-$zmin)*$dimz" | bc -l`
set cscalex = `echo "$dimx/2" | bc -l`
set cscaled = `echo "$dimx*0.8" | bc -l`

#---------------------------------------------------------------
# Use awk to determine the number of grid points along x and z
#---------------------------------------------------------------
set nxz = `awk 'END {print NR-1}' $gridfile`
set nz = `awk 'END {print NR-1}' $profile`
@ nx = $nxz / $nz

####################################
# Velocity vectors
####################################
if ($velocity_vectors == 1) then
  awk 'NR>1 {print "'"$D_scale"'"*$1, "'"$D_scale"'"*$2, "'"$V_scale"'"*$3}' $gridfile >! ux_$tstep
  awk 'NR>1 {print "'"$D_scale"'"*$1, "'"$D_scale"'"*$2, "'"$V_scale"'"*$4}' $gridfile >! uz_$tstep
  set max = `awk 'NR > 1 {if(min==""){min=max="'"$V_scale"'"*$3}; if("'"$V_scale"'"*$3>max) {max="'"$V_scale"'"*$3}} END {printf "%12.5f", max}' $gridfile`
  set vscale = `echo "scale=2; $max/1" | bc -l`
endif

##################
# Temperature 
##################
if ($field == temperature) then
  echo "Plot temperature for time step" $tstep "at time" $ctime
  
  set min = $T_surf
  set max = `echo "$T_surf + $T_scale" | bc -l`

  set xlabel = ""
  set zlabel = ""
  if($dimensional_plot == 1) then 
    set clabel = "Temperature [K]"
  else
    set clabel = "Temperature"
  endif
  set cscale = rainbow
  set ctick = `echo "($max - $min)/4" | bc -l` 
  set stepcontour = `echo "($max - $min)/10" | bc -l`
  set opt = "-Z"
  set adaptive_cscale = 0
  set contour_lines = 1

  set psoutput = temperature_cplot_{$tstep}.ps
  set pdfoutput = temperature_cplot_{$tstep}.pdf
  set jpgoutput = temperature_cplot_{$tstep}.jpg

  awk 'NR>1 {print "'"$D_scale"'"*$1, "'"$D_scale"'"*$2, "'"$T_surf"'"+"'"$T_scale"'"*$6}' $gridfile >! p_$tstep

endif

######################
# Temperature anomaly 
######################
if ($field == tempanomaly) then
  echo "Plot temperature anomaly for time step" $tstep "at time" $ctime
  
  set min = `echo "-($T_scale)/2" | bc -l`
  set max = `echo "($T_scale)/2" | bc -l`

  set xlabel = ""
  set zlabel = ""
  if($dimensional_plot == 1) then 
    set clabel = "Temperature anomaly [K]"
  else
    set clabel = "Temperature anomaly"
  endif
  set cscale = polar
  set ctick = `echo "($max - $min)/2" | bc -l`
  set stepcontour = `echo "($max - $min)/10" | bc -l`
  set opt = "-Z"
  set adaptive_cscale = 0
  set contour_lines = 1

  set psoutput = tempanomaly_cplot_{$tstep}.ps
  set pdfoutput = tempanomaly_cplot_{$tstep}.pdf
  set jpgoutput = tempanomaly_cplot_{$tstep}.jpg

  # store T(x,z) in p1_$tstep
  awk 'NR>1 {print $1, $2, $6}' $gridfile >! p1_$tstep   
  # store <T(z)> in p2_$tstep
  awk -v val=$nx 'NR>1 {Tz[NR]=$3} END {for(k=1;k<=val;k++) {for(i=2; i<=NR; i++) print Tz[i]}}' $profile >! p2_$tstep
  # merge the two files
  paste p1_$tstep p2_$tstep >! p3_$tstep
  awk '{print "'"$D_scale"'"*$1, "'"$D_scale"'"*$2, "'"$T_scale"'" * ($3-$4)}' p3_$tstep >! p_$tstep

endif

######################
# Thermal buoyancy
######################
if ($field == tbuoyancy) then
  echo "Plot thermal buoyancy for time step" $tstep "at time" $ctime
  
  set xlabel = ""
  set zlabel = ""
  set clabel = "Thermal buoyancy"
  set cscale = polar
  set ctick = 0.1
  set stepcontour = 0.05
  set opt = "-Z"
  set adaptive_cscale = 0
  set contour_lines = 1

  set psoutput = tbuoyancy_cplot_{$tstep}.ps
  set pdfoutput = tbuoyancy_cplot_{$tstep}.pdf
  set jpgoutput = tbuoyancy_cplot_{$tstep}.jpg

  # store T(x,z) and alpha(x,z) in p1_$tstep
  awk 'NR>1 {print "'"$D_scale"'"*$1, "'"$D_scale"'"*$2, $6, $11}' $gridfile >! p1_$tstep
  # store <T(z)> in p2_$tstep
  awk -v val=$nx 'NR>1 {Tz[NR]=$3} END {for(k=1;k<=val;k++) {for(i=2; i<=NR; i++) print Tz[i]}}' $profile >! p2_$tstep
  # merge the two files
  paste p1_$tstep p2_$tstep >! p3_$tstep
  awk '{print $1, $2, ($3-$5)*$4}' p3_$tstep >! p_$tstep

  set min = `awk '{if(min==""){min=$3}; if($3< min) {min=$3}} END {printf "%12.5f", min}' p_$tstep`
  set max = `awk '{if(min==""){min=max=$3}; if($3>max) {max=$3}} END {printf "%12.5f", max}' p_$tstep`
  set min = -0.3
  set max = 0.3
  
endif

##################
# Velocity
##################
if ($field == velocity) then
  echo "Plot RMS velocity for time step" $tstep "at time" $ctime

  set min = `awk 'NR > 1 {if(min==""){min="'"$V_scale"'"*$5}; if("'"$V_scale"'"*$5< min) {min="'"$V_scale"'"*$5}} END {printf "%12.5f", min}' $gridfile`
  set max = `awk 'NR > 1 {if(min==""){min=max="'"$V_scale"'"*$5}; if("'"$V_scale"'"*$5>max) {max="'"$V_scale"'"*$5}} END {printf "%12.5f", max}' $gridfile`
  set delta1 = `echo "scale=0; ($max - $min)/10" | bc -l`
  set delta2 = `echo "scale=0; ($max - $min)/5" | bc -l`

  set xlabel = ""
  set zlabel = ""
  if($dimensional_plot == 1) then 
    set clabel = "Velocity [cm/a]"
  else
    set clabel = "Velocity"
  endif
  set cscale = rainbow
  set opt = ""
  set adaptive_cscale = 1
  set contour_lines = 1

  set psoutput = velocity_cplot_{$tstep}.ps
  set pdfoutput = velocity_cplot_{$tstep}.pdf
  set jpgoutput = velocity_cplot_{$tstep}.jpg

  set stepcontour = $delta1
  set ctick = $delta2

  awk 'NR>1 {print "'"$D_scale"'"*$1, "'"$D_scale"'"*$2, "'"$V_scale"'"*$5}' $gridfile >! p_$tstep

endif

##############
# Viscosity 
##############
if ($field == viscosity) then
  echo "Plot viscosity for time step" $tstep "at time" $ctime

  set min = `awk 'NR > 1 {if(min==""){min=log($7*"'"$eta_ref"'")/log(10)}; if(log($7*"'"$eta_ref"'")/log(10)< min) {min=log($7*"'"$eta_ref"'")/log(10)}} END {printf "%12.5f", min}' $gridfile`
  set max = `awk 'NR > 1 {if(min==""){min=max=log($7*"'"$eta_ref"'")/log(10)}; if(log($7*"'"$eta_ref"'")/log(10)>max) {max=log($7*"'"$eta_ref"'")/log(10)}} END {printf "%12.5f", max}' $gridfile`
  set delta1 = `echo "scale=1; ($max - $min)/10" | bc -l`
  set delta2 = `echo "scale=1; ($max - $min)/5" | bc -l`

  set xlabel = ""
  set zlabel = ""
  set clabel = 
  if($dimensional_plot == 1) then 
    set clabel = `echo log@-10@-Viscosity [Pa s]`
  else
    set clabel = `echo log@-10@-Viscosity`
  endif
  set cscale = rainbow
  set opt = "-I"
  set adaptive_cscale = 1
  set contour_lines = 1

  set psoutput = viscosity_cplot_{$tstep}.ps
  set pdfoutput = viscosity_cplot_{$tstep}.pdf
  set jpgoutput = viscosity_cplot_{$tstep}.jpg

  set stepcontour = $delta1
  set ctick = $delta2

  awk 'NR>1 {print "'"$D_scale"'"*$1, "'"$D_scale"'"*$2, log($7*"'"$eta_ref"'")/log(10)}' $gridfile >! p_$tstep

endif

#######################
# Deformation mechanism
#######################
if ($field == defmech) then
  echo "Plot deformation mechanism for time step" $tstep "at time" $ctime

  set min = 0.5
  set max = 5.5

  set xlabel = ""
  set zlabel = ""
  set clabel = "Deformation mechanism"
  set cscale = seis
  set stepcontour = 1
  set ctick = 1
  set opt = "-I"
  set adaptive_cscale = 0
  set contour_lines = 0

  set psoutput = defmech_cplot_{$tstep}.ps
  set pdfoutput = defmech_cplot_{$tstep}.pdf
  set jpgoutput = defmech_cplot_{$tstep}.jpg

  awk 'NR>1 {print "'"$D_scale"'"*$1, "'"$D_scale"'"*$2, $21}' $gridfile >! p_$tstep

endif

###########################
# Strain rate
###########################
if ($field == strain) then
  echo "Plot second invariant of the strain rate for time step" $tstep "at time" $ctime

  set min = `awk 'NR > 1 {if(min==""){min=log($20*"'"$eps_scale"'")/log(10)}; if(log($20*"'"$eps_scale"'")/log(10)< min) {min=log($20*"'"$eps_scale"'")/log(10)}} END {printf "%12.5f", min}' $gridfile`
  set max = `awk 'NR > 1 {if(min==""){min=max=log($20*"'"$eps_scale"'")/log(10)}; if(log($20*"'"$eps_scale"'")/log(10)>max) {max=log($20*"'"$eps_scale"'")/log(10)}} END {printf "%12.5f", max}' $gridfile`
  set delta1 = `echo "scale=1; ($max - $min)/10" | bc -l`
  set delta2 = `echo "scale=1; ($max - $min)/5" | bc -l`

  set xlabel = ""
  set zlabel = ""
  if($dimensional_plot == 1) then 
    set clabel = "log@-10@-(Strain rate) [1/s]"
  else
    set clabel = "log@-10@-(Strain rate)"
  endif
  set cscale = rainbow
  set opt = ""
  set adaptive_cscale = 1
  set contour_lines = 1

  set psoutput = strain_cplot_{$tstep}.ps
  set pdfoutput = strain_cplot_{$tstep}.pdf
  set jpgoutput = strain_cplot_{$tstep}.jpg

  set stepcontour = $delta1
  set ctick = $delta2

  awk 'NR>1 {print "'"$D_scale"'"*$1, "'"$D_scale"'"*$2, log($20*"'"$eps_scale"'")/log(10)}' $gridfile >! p_$tstep

endif

###########################
# Viscous dissipation
###########################
if ($field == viscdiss) then
  echo "Plot viscous dissipation for time step" $tstep "at time" $ctime

  set xlabel = ""
  set zlabel = ""
  set clabel = "log@-10@-\(Viscous dissipation\)"
  set cscale = rainbow
  set opt = ""
  set adaptive_cscale = 1
  set contour_lines = 1

  set psoutput = viscdiss_cplot_{$tstep}.ps
  set pdfoutput = viscdiss_cplot_{$tstep}.pdf
  set jpgoutput = viscdiss_cplot_{$tstep}.jpg

  set min = `awk 'NR > 1 {if(min==""){min=log($8)/log(10)}; if(log($8)/log(10)< min) {min=log($8)/log(10)}} END {printf "%12.5f", min}' $gridfile`
  set max = `awk 'NR > 1 {if(min==""){min=max=log($8)/log(10)}; if(log($8)/log(10)>max) {max=log($8)/log(10)}} END {printf "%12.5f", max}' $gridfile`
  set delta1 = `echo "scale=1; ($max - $min)/10" | bc -l`
  set delta2 = `echo "scale=1; ($max - $min)/5" | bc -l`

  set stepcontour = $delta1
  set ctick = $delta2

  awk 'NR>1 {print "'"$D_scale"'"*$1, "'"$D_scale"'"*$2, log($8)/log(10)}' $gridfile >! p_$tstep

endif

##################
# Composition
##################
if ($field == composition) then
  echo "Plot composition for time step" $tstep "at time" $ctime

  set xlabel = ""
  set zlabel = ""
  set clabel = "Composition"
  set cscale = wysiwyg
  set stepcontour = 0.02
  set ctick = 0.5
  set opt = ""
  set adaptive_cscale = 0
  set contour_lines = 0

  set min = 0
  set max = 1

  set psoutput = composition_cplot_{$tstep}.ps
  set pdfoutput = composition_cplot_{$tstep}.pdf
  set jpgoutput = composition_cplot_{$tstep}.jpg

  awk 'NR>1 {print "'"$D_scale"'"*$1, "'"$D_scale"'"*$2, $9}' $gridfile >! p_$tstep

endif

######################
# Composition anomaly 
######################
if ($field == companomaly) then
  echo "Plot composition anomaly for time step" $tstep "at time" $ctime
  
  set xlabel = ""
  set zlabel = ""
  set clabel = "Composition anomaly"
  set cscale = polar
  set ctick = 0.2
  set stepcontour = 0.2
  set opt = ""
  set adaptive_cscale = 1
  set contour_lines = 0

  set min = -0.4
  set max = 0.4

  set psoutput = companomaly_cplot_{$tstep}.ps
  set pdfoutput = companomaly_cplot_{$tstep}.pdf
  set jpgoutput = companomaly_cplot_{$tstep}.jpg

  awk 'NR>1 {print "'"$D_scale"'"*$1, "'"$D_scale"'"*$2, $9}' $gridfile >! p1_$tstep
  awk -v val=$nx 'NR>1 {Cz[NR]=$6} END {for(k=1;k<=val;k++) {for(i=2; i<=NR; i++) print Cz[i]}}' $profile >! p2_$tstep
  paste p1_$tstep p2_$tstep >! p3_$tstep
  awk '{print $1, $2, $3-$4}' p3_$tstep >! p_$tstep

endif


#########################
# Thermal conductivity
#########################
if ($field == conductivity) then
  echo "Plot thermal conductivity for time step" $tstep "at time" $ctime

  set xlabel = ""
  set zlabel = ""
  if($dimensional_plot == 1) then 
    set clabel = "Thermal conductivity [W/mK]"
  else
    set clabel = "Thermal conductivity"
  endif
  set cscale = rainbow
  set opt = "-Z"
  set adaptive_cscale = 0
  set contour_lines = 1

  set psoutput = conductivity_cplot_{$tstep}.ps
  set pdfoutput = conductivity_cplot_{$tstep}.pdf
  set jpgoutput = conductivity_cplot_{$tstep}.jpg

  set min = `awk 'NR > 1 {if(min==""){min=$12*"'"$k_ref"'"}; if($12*"'"$k_ref"'"< min) {min=$12*"'"$k_ref"'"}} END {printf "%12.5f", min}' $gridfile`
  set max = `awk 'NR > 1 {if(min==""){min=max=$12*"'"$k_ref"'"}; if($12*"'"$k_ref"'">max) {max=$12*"'"$k_ref"'"}} END {printf "%12.5f", max}' $gridfile`
  set delta1 = `echo "scale=1; ($max - $min)/10" | bc -l`
  set delta2 = `echo "scale=1; ($max - $min)/5" | bc -l`

  set stepcontour = $delta1
  set ctick = $delta2

  awk 'NR>1 {print $"'"$D_scale"'"*1, $"'"$D_scale"'"*2, "'"$k_ref"'"*$12}' $gridfile >! p_$tstep

endif

#########################
# Thermal expansivity
#########################
if ($field == expansivity) then
  echo "Plot thermal expansivity for time step" $tstep "at time" $ctime

  set xlabel = ""
  set zlabel = ""
  if($dimensional_plot == 1) then 
    set clabel = "Thermal expansivity [1/K]"
  else
    set clabel = "Thermal expansivity"
  endif
  set cscale = rainbow
  set opt = "-Z"
  set adaptive_cscale = 0
  set contour_lines = 1

  set psoutput = expansivity_cplot_{$tstep}.ps
  set pdfoutput = expansivity_cplot_{$tstep}.pdf
  set jpgoutput = expansivity_cplot_{$tstep}.jpg

  set min = `awk 'NR > 1 {if(min==""){min=$11*"'"$alpha_ref"'"}; if($11*"'"$alpha_ref"'"< min) {min=$11*"'"$alpha_ref"'"}} END {printf "%12.5f", min}' $gridfile`
  set max = `awk 'NR > 1 {if(min==""){min=max=$11*"'"$alpha_ref"'"}; if($11*"'"$alpha_ref"'">max) {max=$11*"'"$alpha_ref"'"}} END {printf "%12.5f", max}' $gridfile`
  set delta1 = `echo "scale=1; ($max - $min)/10" | bc -l`
  set delta2 = `echo "scale=1; ($max - $min)/5" | bc -l`

  set min = 0.35
  set max = 1.85
  set delta1 = 0.2
  set delta2 = 0.2

  set stepcontour = $delta1
  set ctick = $delta2

  awk 'NR>1 {print "'"$D_scale"'"*$1, "'"$D_scale"'"*$2, "'"$alpha_ref"'"*$11}' $gridfile >! p_$tstep

endif

##################
# Heat sources
##################
if ($field == heats) then
  echo "Plot heat sources for time step" $tstep "at time" $ctime

  set xlabel = ""
  set zlabel = ""
  if($dimensional_plot == 1) then 
    set clabel = "Heat production [pW/kg]"
  else
    set clabel = "Ra@-H@-"
  endif
  set cscale = rainbow  
  set opt = ""
  set adaptive_cscale = 1
  set contour_lines = 0

  set psoutput = heats_cplot_{$tstep}.ps
  set pdfoutput = heats_cplot_{$tstep}.pdf
  set jpgoutput = heats_cplot_{$tstep}.jpg
  
  set min = `awk 'NR > 1 {if(min==""){min="'"$heats_scale"'"*$10}; if("'"$heats_scale"'"*$10< min) {min="'"$heats_scale"'"*$10}} END {printf "%12.5f", min}' $gridfile`
  set max = `awk 'NR > 1 {if(min==""){min=max="'"$heats_scale"'"*$10}; if("'"$heats_scale"'"*$10>max) {max="'"$heats_scale"'"*$10}} END {printf "%12.5f", max}' $gridfile`
  set delta1 = `echo "scale=1; ($max - $min)/10" | bc -l`
  set delta2 = `echo "scale=1; ($max - $min)/5" | bc -l`

  set stepcontour = $delta1
  set ctick = $delta2

  awk 'NR>1 {print "'"$D_scale"'"*$1, "'"$D_scale"'"*$2, "'"$heats_scale"'"*$10}' $gridfile >! p_$tstep

endif

#########################
# Depletion
#########################
if ($field == depletion) then
  echo "Plot depletion for time step" $tstep "at time" $ctime

  set xlabel = ""
  set zlabel = ""
  set clabel = "Depletion"
  set cscale = rainbow
  set opt = "-Z"
  set adaptive_cscale = 0
  set contour_lines = 1

  set psoutput = depletion_cplot_{$tstep}.ps
  set pdfoutput = depletion_cplot_{$tstep}.pdf
  set jpgoutput = depletion_cplot_{$tstep}.jpg

  set min = 0
  set max = 1
  set delta1 = `echo "scale=1; ($max - $min)/10" | bc -l`
  set delta2 = `echo "scale=1; ($max - $min)/5" | bc -l`

  set stepcontour = $delta1
  set ctick = $delta2

  awk 'NR>1 {print "'"$D_scale"'"*$1, "'"$D_scale"'"*$2, $19}' $gridfile >! p_$tstep

endif

#########################
# Solidus
#########################
if ($field == solidus) then
  echo "Plot solidus for time step" $tstep "at time" $ctime

  set xlabel = ""
  set zlabel = ""
  if($dimensional_plot == 1) then 
    set clabel = "Solidus [K]"
  else
    set clabel = "Solidus"
  endif
  set cscale = rainbow
  set opt = "-Z"
  set adaptive_cscale = 0
  set contour_lines = 1

  set psoutput = solidus_cplot_{$tstep}.ps
  set pdfoutput = solidus_cplot_{$tstep}.pdf
  set jpgoutput = solidus_cplot_{$tstep}.jpg

  set min = 0
  set max = 1
  set delta1 = `echo "scale=1; ($max - $min)/10" | bc -l`
  set delta2 = `echo "scale=1; ($max - $min)/5" | bc -l`

  set stepcontour = $delta1
  set ctick = $delta2

  awk 'NR>1 {print "'"$D_scale"'"*$1, "'"$D_scale"'"*$2, "'"$T_surf"'"+"'"$T_scale"'"*$18}' $gridfile >! p_$tstep

endif

#######################
# GMT Plot
#######################

#----------------------------------
# Set gmtdefaults
# ---------------------------------
source setgmtdefaults

if($write_time == 1) then
  set title = `echo Time = $ctime`
else 
  set title = ""
endif

#-----------------------------------
# Regridding of the data
#-----------------------------------
surface p_{$tstep} -Gp_{$tstep}.grd -I$xi/$zi -R${xmin}/${xmax}/${zmin}/${zmax}

#----------------
# Make colorbar
#----------------
if ($adaptive_cscale == 0) then
  makecpt -C$cscale -T$min/$max/$stepcontour $opt > p_{$tstep}.cpt
else
  grd2cpt p_{$tstep}.grd -C$cscale $opt -Z > p_{$tstep}.cpt
endif

#--------------------------------------
# Make contour plot and contour lines
#--------------------------------------
grdimage p_{$tstep}.grd -X$X0 -Y$Y0 -R -JX${dimx}/${dimz} -Cp_{$tstep}.cpt -B:${xlabel}:/:${zlabel}:WSen:."${title}": -K >! $psoutput
if ($contour_lines == 1) then
  grdcontour p_{$tstep}.grd -JX -C${stepcontour} -W2/50/50/50 -O -K >> $psoutput
endif

#############################
# Plot velocity vectors
#############################
if ($velocity_vectors == 1) then
  xyz2grd ux_{$tstep} -Gux_{$tstep}.grd -I${deltaux}/${deltaux} -R
  xyz2grd uz_{$tstep} -Guz_{$tstep}.grd -I${deltauz}/${deltauz} -R
  if ($const_vec_leng == 1) then
    grdvector -R ux_{$tstep}.grd uz_{$tstep}.grd -JX -Q0.005/0.1/0.05 -W0.005cm/${vec_col} -G${vec_col} -Sl0.3c -W -O -K >> $psoutput
  else
    grdvector -R ux_{$tstep}.grd uz_{$tstep}.grd -JX -Q0.005/0.1/0.05 -W0.005cm/${vec_col} -G${vec_col} -S{$vscale}c -W -O -K >> $psoutput
  endif
endif

#############################
# Plot streamlines
#############################
if ($streamlines == 1) then
  if (! -f stream) then
    gfortran -o stream stream.f90
  endif
  
  ./stream $dirout $tstep $deltat $nsource $nsteps

  awk '{print $1, $2}' $streamfile | \
  psxy -R -JX -W0.5,${stream_col} -m -O -K >> $psoutput

# Contour plot of the stream function
#  awk '{print $1, $2, $3}' $streamfilet >! p_$tstep
#  surface p_$tstep -Gp_{$tstep}.grd -I$xi/$zi -R${xmin}/${xmax}/${zmin}/${zmax}
#  grdcontour p_{$tstep}.grd -JX -C2 -W1/0/0/0 -O -K >> $psoutput

endif

#############################
# Plot phase boundaries
#############################
if ($phase1 == 1) then
  awk 'NR>1 {print "'"$D_scale"'"*$1, "'"$D_scale"'"*$2, $13}' $gridfile >! p_$tstep
  surface p_$tstep -Gp_{$tstep}.grd -I$xi/$zi -R${xmin}/${xmax}/${zmin}/${zmax}
  grdcontour p_{$tstep}.grd -JX -C0.05 -L0.4/0.45 -W4/${phase_col} -O -K >> $psoutput
endif

if ($phase2 == 1) then
  awk 'NR>1 {print "'"$D_scale"'"*$1, "'"$D_scale"'"*$2, $14}' $gridfile >! p_$tstep
  surface p_$tstep -Gp_{$tstep}.grd -I$xi/$zi -R${xmin}/${xmax}/${zmin}/${zmax}
  grdcontour p_{$tstep}.grd -JX -C0.05 -L0.4/0.45 -W4/${phase_col} -O -K >> $psoutput
endif

if ($phase3 == 1) then
  awk 'NR>1 {print "'"$D_scale"'"*$1, "'"$D_scale"'"*$2, $15}' $gridfile >! p_$tstep
  surface p_$tstep -Gp_{$tstep}.grd -I$xi/$zi -R${xmin}/${xmax}/${zmin}/${zmax}
  grdcontour p_{$tstep}.grd -JX -C0.05 -L0.4/0.45 -W4/${phase_col} -O -K >> $psoutput
endif

if ($phase4 == 1) then
  awk 'NR>1 {print "'"$D_scale"'"*$1, "'"$D_scale"'"*$2, $16}' $gridfile >! p_$tstep
  surface p_$tstep -Gp_{$tstep}.grd -I$xi/$zi -R${xmin}/${xmax}/${zmin}/${zmax}
  grdcontour p_{$tstep}.grd -JX -C0.05 -L0.4/0.45 -W4/${phase_col} -O -K >> $psoutput
endif

if ($phase5 == 1) then
  awk 'NR>1 {print "'"$D_scale"'"*$1, "'"$D_scale"'"*$2, $17}' $gridfile >! p_$tstep
  surface p_$tstep -Gp_{$tstep}.grd -I$xi/$zi -R${xmin}/${xmax}/${zmin}/${zmax}
  grdcontour p_{$tstep}.grd -JX -C0.05 -L0.4/0.45 -W4/${phase_col} -O -K >> $psoutput
endif

#############################
# Plot mesh
#############################
if ($plot_mesh == 1) then
  set xmin = `echo "$xmin / $D_scale" | bc -l`
  set xmax = `echo "$xmax / $D_scale" | bc -l`
  set zmin = `echo "$zmin / $D_scale" | bc -l`
  set zmax = `echo "$zmax / $D_scale" | bc -l`
  psxy -JX -R${xmin}/${xmax}/${zmin}/${zmax} -W1,0/0/0 -m $meshfile -O -K >> $psoutput
endif

#############################
# Plot colorbar
#############################
psscale -D$cscalex/-0.5/$cscaled/0.5h -Cp_{$tstep}.cpt -N -Ba${ctick}:"${clabel}": -O >> $psoutput

####################################
# Conversions, moving and cleaning
####################################

ps2raster $psoutput -Au -P -Tf -F$pdfoutput
if ($makejpg == 1) then
  ps2raster $psoutput -Au -P -Qt4 -Qg4 -E$jpgdensity -Tj -F$jpgoutput
  mv -f $jpgoutput $dirplots$field
endif

if ($makejpg == 1 && $deletepdf == 1) then
  rm -f $pdfoutput
else 
  mv -f $pdfoutput $dirplots$field
endif

if ($streamlines == 1) then
  mv -f $streamfile $dirout
endif

rm -f p*_{$tstep} $psoutput p*_{$tstep}.grd p*_{$tstep}.cpt ux_{$tstep}* uz_{$tstep}*


